Elementary Excitation Modes in a Granular Glass above Jamming 
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The dynamics of granular media in the jammed, glassy region is described in terms of "modes", 
by applying a Principal Component Analysis (PCA) to the covariance matrix of the position of 
individual grains. We first demonstrate that this description is justified and gives sensible results in 
a regime of time/densities such that a metastable state can be observed on long enough timescale to 
define the reference configuration. For small enough times/system sizes, or at high enough packing 
fractions, the spectral properties of the covariance matrix reveals large, collective fluctuation modes 
that cannot be explained by a Random Matrix benchmark where these correlations are discarded. 
We then present a first attempt to find a link between the softest modes of the covariance matrix 
during a certain "quiet" time interval and the spatial structure of the rearrangement event that 
ends this quiet period. The motion during these cracks is indeed well explained by the soft modes 
of the dynamics before the crack, but the number of cracks preceded by a "quiet" period strongly 
reduces when the system unjams, questioning the relevance of a description in terms of modes close 
to the jamming transition, at least for frictional grains. 

PACS numbers: 



I. INTRODUCTION 

The physical processes by which super-cooled liquids, 
granular systems, and colloids acquire rigidity are not 
well understood. At first sight, the phenomenon of rigid- 
ity is utterly trivial: we know that when we move one 
end of a ruler the other end moves the same distance. It 
is so simple that, as P.W. Anderson noticed, it is hard to 
realize that such an action at a distance is not built into 
the laws of nature except in the case of the long-range 
forces such as gravity and electrostatics... We are so ac- 
customed to this rigidity property that we don't accept its 
almost miraculous nature, that is an "emergent property" 
not contained in the simple law of physics, although it is 
a consequence of them. 

Recently, intense research has been devoted to this 
problem and it has become clear that the emergence of 
rigidity in soft matter is likely to be related to a collec- 
tive phenomenon. Many hints came from numerical and 
analytical studies of the jamming transition of hard and 
elastic frictionless spheres [ll, Q . In this case, it has been 
shown that when the system acquires rigidity it has no 
redundant mechanical constraints. As a consequence, it 
is in a marginally stable, isostatic, state. This has dra- 
matic consequences for the vibrational spectrum, which 
displays a broad band of soft modes [1, Ql • The role of 
these modes in the dynamics close to the rigidity tran- 
sition and, more generally, for glassy liquids has been 
emphasized in [sHll- However, the applications and ver- 
ifications of these theoretical ideas in experiments are 
scarce. A first attempt has been performed for colloidal 
glasses Q, but the limitation in experimental resolution 
does not allow to draw definitive conclusions. Here we 



focus on mechanically driven granular media. These are 
the physical systems that triggered the studies of anoma- 
lous properties of vibrational modes and isostatic prop- 
erties 0^ [l^. Despite of this, there is still no experi- 
mental study in the literature on the role of the modes 
close to the rigidity transition. The aim of our work is to 
present a first analysis of the modes close to the rigidity 
transition of vibrated frictional grain assemblies. Note 
that the presence of friction is expected to modify the 
properties of the transition compared to the ideal case of 
hard spheres [ill, [l^. In particular, our system seems to 
be characterized by micro-cracks of all scales, leading to 
'jumps' in the position of particles with a power-law dis- 
tribution of sizes [l^ , which makes the analysis in terms 
of modes particularly tricky. Still, we believe that the 
tools we developed are interesting also from a method- 
ological point of view, and will be useful for analysing 
other systems that undergo a jamming transition. 



In [T3 | it was shown that as the packing fraction of a 
horizontally vibrated monolayer of bidisperse hard grains 
is increased beyond a certain packing fraction (f>j, the sys- 
tem is able to support mechanical stresses. This is the 
rigidity transition, which appears as a genuine critical 
point, where a dynamical correlation length and a corre- 
lation time simultaneously diverge, showing that the dy- 
namics occurs by involving progressively more collective 
rearrangements 4>j. Contrary to the case of frictionless 
hard sphere or colloids the pressure does not diverge at 
but at a higher density. 



Experimentally, we have access to the covariance ma- 
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trix of the positions, Cp ^. Whether and to what extent 
this can be interpreted in terms of vibrations along some 
modes is one of the main open questions that we shah 
address. We shah also investigate how the eigenstates 
and eigenvalues of Cp evolve when approaching and 
their relation with the dynamics. In order to do that, we 
have to separate signal from noise in the cigenproperties 
of Cp. This is a common and crucial problem in dealing 
with covariance matrices, which will be addressed by us- 
ing tools and concepts previously developed and used in 
other fields like finance and biology |15l - [20| . 



II. EXPERIMENTAL SYSTEM AND 
PRELIMINARIES ON THE PARTICLE 
POSITIONS COVARIANCE MATRIX Cp 

The experimental set-up and the quench protocols are 
described in detail in [ij]. A 1:1 bidisperse monolayer 
of 8500 brass cylinders of diameters ds = 4 ± O.Olmm 
and db — 5 ± O.Qlmni stands on a glass plate which is 
horizontally vibrated at a frequency of 10 Hz and an am- 
plitude of lOmm. The grains are confined within a fixed 
rectangular metal frame of width L ?» 100 ds- The pack- 
ing fraction (/) can be adjusted by moving a lateral wall 
on which we control the pressure. The stroboscopic mo- 
tion of a set of 1500 grains in the center of the sample 
is tracked with an accuracy of 2.10~^ds. Lengths are 
measured in ds units and time in cycle units. The ini- 
tial protocol produces a very dense state with a packing 
fraction of (j) = 0.8457. The packing fraction is then de- 
creased by very small steps down to 0.84. For each (j), 
the plate vibrates 10^ cycles during which the pressure 
at the wall is stored. At high packing fraction, the mean 
pressure is dominated by the static pressure, which is 
measured by interrupting the vibration. At some 4>, the 
kinetic part of the pressure becomes dominant and this 
is identified as the jamming transition, which takes place 
at e [0.8417,0.8422]. 

The main properties of the grain displacements have 
been discussed in detail in [l3, [2l| . Very recently, we re- 
examined these statistics of the displacements and found 
the rather surprising results alluded to above, which we 
report in another paper of the present special issue p^ . 
First, let us insist on the fact that the typical displace- 
ment of the particles is of the order of one hundredth of 
its diameter. Accordingly, all structural rearrangements 
are frozen on the experimental timescales: the neighbors 
of a given particle do not change during the experimen- 
tal time scale. Second, the motion of the particles, sub- 
diffusive at short times and diffusive at asymptotically 



We have also studied the covariance matrix of the instantaneous 
velocities but at the present stage of the study, it did not provide 
further insight. We thus concentrate here on the results given by 
the study of Cp 



large times, exhibits a super-diffusive motion at interme- 
diate timescales close to the jamming transition. Our 
recent analysis of the data shows that this superdiffu- 
sion is not induced by long-range temporal correlations 
of the velocity field, as we first surmised in [13. Quite 
on the contrary, the displacements on the intermediate 
timescale are made of a large number of incoherent jumps 
with a broad distribution of jump sizes. However, these 
jumps become more and more collective as the systems 
becomes rigid at (pj, which appears as a genuine critical 
point, where the dynamical correlation length diverges. 

As stated in the introduction, our aim here is to further 
characterize the dynamics and its spatial organization 
close to the rigidity transition by studying the covariance 
matrix of the particles positions, defined as: 

Cp = ( Sri,a drj,j3 )t = ( (n,Q-(rj^Q)T) {rj,f3~{rj,p)T) )t, 

where r.i^a is the a = x or y Cartesian coordinate of the 
^th gj-j^jjj a^jjfj { ■ )t denotes the temporal average over an 
observation window of duration T. 

For solids at thermal equilibrium, the modes of Cp can 
be identified with structural vibrational modes because 
particles simply oscillate around their equilibrium posi- 
tions. For example for crystals at low enough tempera- 
ture the matrix Cp is equal to the temperature times the 
inverse of the Hessian matrix of the potential energy eval- 
uated for the ground state configuration. In this case the 
eigenvectors of Cp are plane waves that identify with the 
phonons. In the present case, the system being driven out 
of equilibrium, it is not warranted at all that the modes 
of Cp can be interpreted as vibrational modes. However, 
as argued in the introduction and as will be confirmed in 
the following, studying the spectral properties of Cp re- 
mains a powerful tool of investigation, provided that the 
particles have a well defined average position: Cp mea- 
sures the fluctuations around a metastable state and its 
spectral properties allow one to interpret these fluctua- 
tions in terms of effective excitation modes. 

We thus start by investigating the fluctuations of the 
particle positions around the average position. Here we 
focus on a single component x of the position, but we 
have checked that the conclusions are identical for both, 
confirming that the dynamics is isotropic as already ob- 
served in \XM- -f'or a given particle i, one can compute the 
average position {xi)T on a time T and the fiuctuations 
around it Sxi(t) = Xi{t) — {xi)T- Note that in glassy dis- 
ordered systems, this average position can itself evolve 
with time and be an extra source of fluctuations. The 
variance of 5xi over time T characterizes how far the par- 
ticle is, typically, from its average position: af = {Sxf)T- 
We find (see below) that ct; significantly fluctuates from 
particle to particle, reflecting the presence of dynami- 
cal heterogeneities in the system: while some particles 
hardly move during time T, others are able to "rattle" 
quite a bit (but still on scales much smaller than the 
grain diameter!). More precisely, the distributions p{(Ti) 
arc shown on the right of Fig. [TJ When decreasing the 
packing fraction towards p{ai) shifts to larger values 
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FIG. 1: Left: Distributions of tlie position fluctuations p{5xi/ai) for four values of the packing fraction and four durations of 
observation T. We focus on the top of the distribution, which is compared to a Gaussian (continuous line). Right: (a) p{(Ji) 
for different cj> keeping T = 100 constant. The arrow indicates the direction of decreasing (f) e [0.8417, 0.8426, 0.8440, 0.8457] 
(b) p((Ti) for different T keeping </> = 0.8426 constant. The arrow indicates the direction of increasing T e [10^10^10*]. 



of CTi, indicating larger overall motions for eaeh particle, 
as expected. As (f> decreases, p((Ti) also broadens sig- 
nificantly demonstrating more and more heterogeneities 
among the particles. Indeed describing the right tail of 
the distribution by a power law: p{(7i) ~ ct~^~^, one 
find p, decreasing from k, 4 to ~ 3, when decreasing 
towards 0j. As a matter of fact, for the largest pack- 
ing fractions, the power-law tail is so steep that it can 
as accurately be described by an exponential. When T 
increases, the distribution p{ai) shifts to larger values of 
(Ti as expected, but does not broaden, indicating that the 
heterogeneities are already well developed within an in- 
terval of time T = 100. Such observations are yet another 
confirmation of the statistical properties of the dynam- 
ics studied in [ij, [13] • Note that the exponent /i here 
should not be confused with the exponent describing the 
tail of individual jump sizes, as defined in [l^: here, we 
characterize the variation of the vibrations across differ- 
ent grains, and not for a single grain over time. In order 
to perceive the difference more clearly, imagine a case 
where all particles perform exactly the same motion, be 



it a regular random walk or a Levy flight: in both cases, 
p{(Ti) should then be a delta function since there is no 
dispersion at all. 

We then study the distribution of rescaled positions, 
5xi/ai, by averaging over all times and all particles. The 
distributions are computed for four different packing frac- 
tions (p s-nd four durations T of the window of observa- 
tion. They are then ensemble averaged over the lO^'/T 
intervals provided by the full dataset. From now on all 
statistical quantities (such as the eigenvalue spectra, etc.) 
are evaluated this way, without further specifying it ex- 
cept when necessary to avoid confusion. 

The distributions shown in Fig. [T] highlight some im- 
portant characteristics of the dynamics. The parameter 
space ((/), T) can be divided into two regions, as illus- 
trated by the hatched line: for small enough observation 
duration T or large enough packing fractions, the dis- 
tributions are unimodal with a Gaussian core: particles 
jiggle around a well defined average position; for longer 
T or smaller densities, the distribution starts developing 
a flat top, with a poorly deflned maximum. This sug- 
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gest that on these longer observation times, the average 
position of a significant part of the particles is not well de- 
fined anymore. Particles either drift slowly or even find 
(collectively) another metastable position, as suggested 
by the double peak observed in the case 4> = 0.8417 
and T = 10^, i.e. for the loosest packing fraction and 
the longest observation time. This means that over long 
time scales, the evolution of the average position becomes 
comparable or even larger than the fluctuations, and it 
becomes meaningless to describe the system in terms of 
small vibrations around a fixed metastable state. For an 
infinite size system, some rearrangement always happens 
somewhere, and the covariance matrix Cp is always ill- 
defined. The "allowed" time scale Tniax{4>, L) is expected 
to scale inversely with the system size; however, when 
Tmax becomes too small, statistical noise becomes domi- 
nant and prevents a reliable estimation of the spectrum 
of Cp. In the following sections, wc will navigate be- 
tween these constraints and try to identify well defined 
eigenmodes of the motion. 
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FIG. 2: Spectral properties for the whole system {<j) — 0.844, 
N — 1550 particles, T = 100). (a): Normalized spectrum 
Am VS. m for the experimental data and for the Random 
Matrix case (RMcr)- The (blue) dotted line with exponent 
— 1/a = —1 is the prediction for the 2D Crystal, (b): The 
associated density of states. Dashed lines are guide for the 
eyes. 



A. The role of noise 



III. SPECTRAL PROPERTIES OF Cp 



In this section we shall study in detail the spectral 
properties of Cp. Our aim is twofold: first, as stated 
above, the spectral properties are affected by measure- 
ment noise for finite T. Thus it is important to disen- 
tangle trivial properties of Cp induced by the noise from 
relevant ones, which we do in the first subsection. Sec- 
ond, we would like to understand whether Cp is indeed 
measuring some steady fiuctuations around a well-defined 
metastable state. We will refer to such a property as ro- 
bustness and study it in the next subsection. Third, we 
will see that the structure of the modes itself confirms 
that the 10 first modes are significantly out of the noise 
range. In order to do so we concentrate on one spe- 
cific case in the middle range of our parameter space, 
4> = 0.844 and T = 100, for which particle positions 
seem to be well defined on the observation window du- 
ration and we consider the whole set of tracked particles 
N = 1500. Given that the correlation matrix is com- 
puted in a observation window T < 2N, there arc at 
the best only T non-zero eigenvalues amongst 2N. The 
eignevalues are normalized by (t^/Q, where a — {ai)i is 
the average of the di's over all particles and Q = T/2N 
is the total number of measured data points 2N x T di- 
vided by the total number of variables AN'^ . With such a 
normalization, one can easily compare the spectra of Cp 
for systems with different average mobility a as well as 
for computations of Cp with different values of Q ~ for 
instance when considering subsystems of smaller size A^, 
as we shall do in the next section. 



In order to obtain some hints on the role of noise in the 
spectral properties of Cp we will compare our results to 
the ones obtained by constructing the covariance matrix 
with iid random variables Zi{t) = r]i{t)ai, where rjiit) are 
iid Gaussian variables and at are positive random vari- 
ables following the experimental distribution p{(Ji). Note 
that cr, is constant during each interval of duration T . 
In this benchmark model, which we shall refer to as the 
Random Matrix case {RMa), the spatial correlations be- 
tween rii{t) and r]j(t), j ^ i, and between the different 
cr's are discarded. By comparison, we will be able to 
evaluate the relevance of these correlations in the experi- 
mental system. Our results will also be compared to the 
case of a 2-d equilibrium crystal. 

Figure [2] displays the normalized spectrum Am and the 
associated density of states p{X) for both the experimen- 
tal data and the RMa- simulation when (j> = 0.844 and 
T = 100. Note that there is a straightforward corre- 
spondence between the behaviour of the more commonly 
studied p(A) at large A and that of Am at small m since 
-^m(A) is exactly the inverse cumulated distribution of 
A. Accordingly a powerlaw behaviour Am. m~^/" trans- 
lates into a density p{X) ^ A~'^+"-'. For instance a 2D 
crystal, with a density of states p{uj) ^ w'^"^ = lu, with 
u} ^ 1/A^/^, has p(A) ~ A~^, that is a = 1 and Am "^"^ 
as indicated on the figure by the blue dotted line. For the 
Random Matrix case (RMa-), if the distribution p(cri) has 
power-law tails with exponent p., then the top eigenvalues 
of the correlation matrix also has a power-law tailed dis- 
tribution, with exponent l-t-^/2 and p(A) should decay at 
large A at least as slow as A^'-^^^^^'. In the present case, 
/i « 4 and one would expect a ~ 2, whereas we mea- 
sure here 1/a ~ 1/4. The reason for this discrepancy 
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is insufficient sampling: as is clear from Fig. [1] right, 
there is hardly a factor 10 contrast between the largest u 
encountered in the sample and the typical one as given 
by the median of the distribution. We have check that 
with many more samples, the expected power-law tail 
eventually appears. But we have been careful here to 
take exactly the same statistics in the simulation and in 
the experiment, so that the comparison made in Fig. [5] 
is meaningful. The ten largest eigenvalues for the ex- 
perimental system are therefore clearly larger than for 
the Random Matrix case. Our data is consistent with 
a ~ 2/3, that is a slower decay of the spectrum than 
for both the RM^, and the crystal case. This compari- 
son shows unambiguously that the top eigenvalues of Cp 
contain useful information about the dynamics of the sys- 
tem, and arc not drowned in noise. It also demonstrates 
the existence of strong spatial correlations: by moving to- 
gether, particles achieve large collective fluctuations that 
would not develop otherwise. 

B. Micro-cracks and robustness 

As already stressed, in such an heterogeneous system 
the above analysis strongly relies on the selection of the 
observation windows, in order to ensure that the system 
remains in a single metastable state. We thus compute 
the instantaneous self density-correlation function: 

C^{tM) = (cos(q'- -r,(to)]))», (1) 

where (t) is the particles position at time t and (f is a 
wave vector whose amplitude is given by g = 7r/a. a is 
chosen as a small length scale of the order of a* = ((f^ {t+ 
T*) - ri[t)Y)^/'^ = 7. 10~^, where r* is the timescale at 
which dynamical heterogeneities are maximal (see [l3j in 
the present volume for details). The average is computed 
over all particles, but not on the initial time and it is 
not ensemble averaged either. Cq(t,tQ) decays to zero 
when the average displacement of the particles between 
time to and to + t is larger than a. 

Ideally, one would like to observe sudden drops of 
Cq(t, 0) that signal moments when a significant collective 
event occurs, hopefully separated by long enough "quiet" 
periods. Also, the same plateaus and cracks should be 
present for a reasonable range of length-scales a. This is 
not the case here, as clearly observed in Fig. EJa). For 
all a G [a*, 10a*], the decrease of Cq{t,0) is progressive 
rather than taking place during sudden drops. Further- 
more, the observation time windows which look like quiet 
for a given a, are in fact very jerky when decreasing a. 
The reason for these features are (i) the heterogeneity 
of the relaxation (for a large enough system, some re- 
laxation event is always taking place somewhere in the 
system) and (ii) the scale invariance of these relaxation 
events, or "micro-cracks" as recently pointed out in [l3| . 

This scale invariance makes it very hard (if not im- 
possible) to define properly the metastable states of the 
system, and the corresponding covariance matrix Cp. It 
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FIG. 3; Relaxation events for (j> = 0.844. Top: (a): Instanta- 
neous self density-correlation function Cq{t,to) as defined in 
the text. Each curve is for a different a £ [a*, 2a*, 5a*, 10a*] 
and the arrow indicates increasing a. (b) positions of the 
particles which where identified as significantly contributing 
to the relaxation of the system (see definition in the text). 
Different symbols (and colors online) indicate different time 
windows of duration T = 1000. The square box surrounds 
the rather inactive region, which we select as the subset of 
particles over which we re-compute Cq{t,to) displaid in (c). 
The dotted lines indicate the temporal windows isolated as 
quiet periods, (d): Same plot as in (b) but restricted to the 
period of time t G [3000 6000]. 

suggests to identify relaxation events not by their size 
but through an iterative process such as the one pro- 
posed in [2^ to identify cage jumps in the trajectories 
of particles in a super-cooled liquid. In a nutshell the 
algorithm consists in cutting each trajectory in two sub- 
trajectories, in such a way that each subset maximizes a 
clustering criteria, and in applying iteratively the algo- 
rithm to each subset until the maximization criterion is 
no more significant. As a result one obtains for each par- 
ticle a set of instants corresponding to the times when it 
has relaxed, without specifying the amplitude required to 
relax. Fig. ^h) displays such events: at all times some 
small regions of the system relax, each of them contribut- 
ing to a small decrease of Cq{t,0). 

Altogether, when approaching the jamming transition 
from above, the system as a whole becomes more and 
more heterogeneous, less rigid, and metastable states 
harder and harder to define. This makes the computation 
of Cp increasingly difficult, precisely where we would like 
to use it to characterize the dynamics. However, one also 
notices in Fig. [3jb) a region indicated by the square box, 
where there is little activity as compared to the rest of 
the system. Fig.|31Jc) again displays Cq{t, 0) but averaged 
on the particles belonging to this quiet region only. One 
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FIG. 4: Robustness of the modes for — 0.844 as charac- 
terized by R{m) for T = 100 and T = 1000; (a): whole set 
of particles, (b) subset of particles and observation window 
identified in fig.[3](iV = 350 particles, T = 100 and T = 1000). 



now can better identify periods of time where Cg(t, 0) is 
rather constant, independently of a. These sub-regions 
are rigid during long enough time intervals to perform 
the analysis in terms of modes. In the following we shall 
refer to these sub-systems and time-interval as the "rigid 
subsets" of the system. 

This we confirm by assessing the robustness of the 
modes. For that purpose, we compute the following in- 
dicator: 



3=+M 

R{m)= {K 



\^'m+j) 



(2) 



where (Am|A'^^^) is the scalar product between the 
modes computed during two successive observation win- 
dows of duration T. If the two eigenbases are precisely 
the same, R{m) is equal to 1 for all eigenvectors A,„. Note 
that the definition allows that the modes computed in 
one observation window project onto any of the (2M-I-1) 
modes of the second basis surrounding mode m, in order 
to allow the neighbouring modes to possibly exchange 
their rank. For M > 2 and not too large the results 
are basically independent of M. Here we fixed M = 2. 
One observes in Fig. |3] that the robustness of the modes 
is twice larger when restricting the analysis to the rigid 
subsets. 

Let us finally describe the spectral properties within 
the rigid subsets. Figure© displays the distribution 
p{c7i) and the spectrum Am, which we compare to the 
ones obtained for the whole system (figHJ-right and[2]-a). 
We shall come back to the description of the distributions 
p{X) in the next section. The right tail of the distribu- 
tion p{<7i) is more narrow (p ~ 6) for the subset than 
for the whole system (/x ~ 4) confirming that these rigid 
subsets are more homogeneous. For the RM^ case, the 
spectrum A™ for small m is again fiatter than expected, 
due to insufficient sampling (l/a ~ 1/4 instead of 1/3). 



FIG. 5: Spectral properties computed for the subset of par- 
ticles and observation window identified in fig. [3] ((^ = 0.844, 
iV = 350 particles, T = 100). (a): Density of ai for both 
the whole system and the subset of interest (b): Normalized 
spectrum Am vs. m for the experimental data and for the 
Random Matrix case [RMc,)- The (blue) dotted line with 
exponent —1/ot = —1 is the prediction for the 2D Crystal. 
Dashed lines are guide for the eyes. 



More remarkable is the fact that the spectrum for the 
experimental system remains well above the RM„ case 
and that it is almost identical to the one obtained in the 
whole system, suggesting that (i) it is not dominated by 
the shape of the distribution p{(Ji) but on the contrary 
unveils non trivial correlations; (ii) the heterogeneities 
associated to these correlations are present at all scales. 

Altogether, despite rather poor statistics and a signif- 
icant amount of noise in the spectral properties of Cp, 
the difference reported between the random matrix and 
the experimental cases confirms that one can trust the 
largest eigenvalues and that the spectrum in the experi- 
mental system is mostly governed (in its top region) by 
non-trivial spatial correlations. 



C. The structure of top eigenvectors 

We can substantiate this last assertion even more by 
comparing the localization properties of the associated 
eigenvectors. We start with the RM„ case. Assum- 
ing that the Ui are power law distributed with an ex- 
ponent p, the maximum Umax is given by the equation 
N Jq""" p{ai)d(7i ~ 0(1), which leads to (J,nax k N^/f^. 
Calling i* the value of i corresponding to the maximum 
(T;, one expects in the absence of correlations that the 
covariance matrix has a very large diagonal entry at i* 
(in the present case, the two eignevalues corresponding 
to the X and y directions are equally large as imposed 
by the equally large ai in both directions). A reason- 
able guess, that can be justified using arguments as the 
one developed in [2^, is that this leads to the largest 
eigenvalue and that the rest of the covariance matrix can 
be considered as a perturbation. Accordingly the largest 
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FIG. 6: Participation ratio of the modes computed within a 
rigid subset for both the experimental data and the Random 
Matrix case (RM^); cj> = 0.844. 



eigenvalue in the RM^ is given by cr^j and the corre- 
sponding eigenvector is completely localized on i* . We 
have checked that both facts are indeed very well real- 
ized. Note that similar results hold for the second, third,., 
largest eigenvalues which are related to the second, third 
largest value of ai . 

When considering the experimental covariance matrix, 
the largest eigenvalues and their corresponding eigenvec- 
tors are instead very different. This allows us to make 
clear that they are not due at all to noise and that spa- 
tial correlations are very instrumental in creating large 
eigenvalues or large fluctuations in the particle positions. 
In order to quantify this effect, we compute the partici- 
pation ratio defined as: 



Pirn) 



1 



(3) 



where -itj^ is the normalized displacement of particle i 
within mode m. This quantity is such that, if the mode 
m is completely localised on one particle, P{m) = 1/N. 
The other extreme case is when all the particles con- 
tribute equally to the mode: in this case P{m) = 1. 
Fig. ini displays P{m) for both the experimental and the 
Random Matrix cases. It is clear that largest eigenvalues 
for the RM„ case have a very small participation ratio, as 
expected since they are essentially localized on one or a 
few sites. Instead, the modes from experiments are char- 
acterized by a much higher P{m), indicating that these 
modes are delocalized although less than plane waves for 
which P{m) = 2/3. Beyond m = 10, the participation 
ratio in both cases are very similar, suggesting that these 
bulk modes are incoherent and dominated by the local 
fluctuations of Ci, and not by spatial correlations. 

After this long but necessary description of the 
methodology, we can venture into the investigation of 
the relevant eigenvalues and modes structure, when ap- 



proaching the jamming transition. On the basis of the 
above analysis, we will concentrate the computation of 
Cp on the rigid subsets and restrict the analysis to, say, 
the 10 largest eigenvalues and corresponding eigenvec- 
tors. 



IV. TOWARDS THE JAMMING TRANSITION : 
MODE STRUCTURE AND DYNAMICS 

In the following, we present a quantitative analysis of 
the modes and their properties when approaching the 
rigidity transition. This is to our knowledge the first 
attempt of this kind for granular assemblies. We first 
characterize the mode structure and then assess the role 
of these modes in the dynamical evolution of the system. 

We focus on the spectral properties of Cp for the four 
densities G [0.8417, 0.8426, 0.8440, 0.8457] following the 
methodology outlined in the previous section, i.e. iden- 
tifying rigid subsets for which we measure Cp. Then 
for a given density, we average over all the available 
rigid subsets. For the two densest cases, we identified 
two sub-regions which are rigid during typically 3000 cy- 
cles. Closer to the jamming transition, there is no re- 
gion, which remains rigid during more than 400 cycles. 
We identified 6 of such rigid subsets for = 0.8426 and 
5 for = 0.8417. In all cases the regions have about 
N = 350 particles. An important observation is that it 
becomes increasingly difficult to measure the modes when 
approaching the rigidity transition. This is likely related 
to the findings explained in the companion paper [isj 
which show that at the dynamics is due to temporally 
incoherent but spatially correlated Levy jumps, corre- 
sponding to micro-cracks of all amplitudes that span the 
system, making it hard to find sub-regions where nicely 
separated, "big" cracks occur. 



A. Structure of the modes close to (pj 

The study of the spectral properties when approaching 
unveils that the softest modes become both softer and 
more extended: 

• Mode softness: As observed on figure [71 the largest 
eigenvalue increases when approaching more- 
over Am vs. m becomes steeper in the log-log plot 
leading to an exponent 1 /a slowly varying between 
3/2 and 2. Accordingly the spectrum develops 
larger tails and there is a redistribution of spectral 
weight towards larger eigenvalues. The situation is 
clearly different from the crystal case were 1/a = 1. 

• Mode extension: The closest the system is to the 
jamming transition, the more coherent and spa- 
tially organized the softest mode is. This is visually 
clear on Fig. [51 which provides an example of the 
softest mode for each packing fraction. 
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FIG. 7: (a) Normalized spectrum Am vs. xm for (f) £ 
[0.8457,0.8440,0.8426,0.8417]. As specified in the text, Cp 
is computed within the subset of particles (A'' ~ 350), with 
well identified rigid periods, (b) Corresponding eigenvalue 
densities p(Am). 



Two quantitative results support this last assertion. 
First, the participation ratio for the largest modes in- 
creases from P(l) = 0.2 to P(l) = 0.4 when approach- 
ing (pj. Second, spatial correlations within the modes 
increase. This is measured by computing the following 
correlation function: 



Cm{r) = - ('a™)).K„ - (Um))) 



(4) 



where the average is computed over all pairs of particles 
separated by r. These spatial correlators are plotted for 
m = 1 and m = 10 in the insets of Fig. [H Clearly, 
the correlation extends on a longer distance when the 
system is closer to cfij. An interesting feature is that 
Cmir) becomes negative for r « 10, indicating some anti- 
correlation, which we attribute to the vortices pattern 
observed in the modes. Also, the correlation is much 
weaker for m = 10 than for m = 1. This effect is further 
characterized in the main plot of the same figure, where 
we plot Im = J2r<5 ^m('') versus m for the four packing 
fractions. Not only the modes have a structure on a larger 
scale closer to (f>j, but also more of them are structured. 



B. Soft modes and dynamics in a metastable state 

We now turn to the relation between the modes |A„i) 
and the dynamics \r{t)) = {ri{t)}. We first concen- 
trate on the dynamics restricted to the rigid subsets 
and consider the projection of the real dynamics on 
the modes computed in an observation window preced- 
ing the dynamics by a lag time r. More precisely, let 
[t, t + T] be the time window where the basis of eigen- 
modes {|Am)} is computed. The dynamical evolution 
\Sr{T)) — \r{t + T + r)} — \r{t + r)) is then projected 
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FIG. 8: Top: Four realisations of the first mode (m = 1) 
for each packing fraction (ji G [0.8457,0.8440,0.8426,0.8417]. 
Bottom: Spatial correlations within the modes. The main 
plot provides an estimation of the correlation length as a func- 
tion of the mode rank for the four packing fractions. Insets: 
spatial correlation for the four packing fractions for the modes 
m = 1 and m = 10. 



on the modes m and the corresponding component is 
rescaled by the amplitude of the dynamics: 



(A,„|,5r(T)} 

Cm — 



(fr(T)|5r(T)) 



(5) 



The components Cm satisfy X]m('^™)^ ~ 1 since the eigen- 
vectors form a complete basis. We sort the Cm in decreas- 
ing order c^. = ci > C2... > and, following [Bj, define: 



F{m) 



k=l 



(6) 



F{m) measures the fraction of the dynamics "explained" 
by the m most contributing modes. Here, in the light of 
the previous section, we have chosen to consider the 10 
first modes, T = 100, t varies from 1 to 1,000 cycles and 
we average J^(IO) on two initial times t as well as over 
all the rigid subsets. Figure O-top displays F{1Q) for the 
four packing fractions. Three key aspects emerge: 
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• For all packing fractions F{10) fluctuates around a 
large constant value. This shows that even for large 
T the dynamics is well described by the 10 most sig- 
nificant modes of Cp as long as the system remains 
in a metastable state. From this perspective, the 
modes defined by Cp for the rigid subsets give a 
faithful representation of the dynamics and can be 
indeed considered as effective vibrational modes. 

• Interestingly, the average value {F{10))r increases, 
beyond error bars, when (j) increases towards (j)J- 
This indicates that the 10 first modes concentrate a 
more important part of the dynamics as (/> — > 0j, in 
agreement with the idea that the dynamics becomes 
more collective, or structured. This is further 
demonstrated on Figure [3]-bottom where {F{m,)) 
is plotted versus m for the fifty largest modes. The 
closer to (/>,/, the larger is {F{m)). 

• The fiuctuations of F{10) are clearly more corre- 
lated in time when (j) decreases, revealing that the 
modes have a larger characteristic time closer to 



C. Soft modes and cracks: a preliminary analysis 

The observations above suggest that when approach- 
ing a smaller and smaller amount of modes concen- 
trate the dynamics for longer and longer times : the sys- 
tem "rattles" around its metastable state along more and 
more preferential and softer directions in phase space. 
However, for all the methodology issues seen above, and 
because of the lack of theoretical grounds for frictional 
systems, determining whether such directions determine 
the way the system locally cracks (as seen in other sys- 
tems 0, [131 ) is an extremely challenging issue in the 
present system. Here we provide a first attempt to answer 
this question, which obviously deserves further analysis. 

We identify in the self-density correlation function 
Cq{t) a sudden drop preceded by a quiet period (see 
Fig. rroi-a). We compute the covariance matrix of the 
positions restricted to this rigid subset before the crack 
r< ~ [ti,t2]. We then consider two further time intervals 
related to the crack: one during the crack T= = [^2,^3] 
and one after the crack T> = [t3,t4]. We observe the 
real dynamics during these intervals and project the dy- 
namics onto the first m modes determined before the 
crack, defining (m) , _FU (m) and f > (m) . One observes 
on figure [TU]-(b) that F^{m) and F^{m) share a simi- 
lar behaviour as a function of m, different from the one 
of Fy{m): both F^{m) and F^{m) increase sharply at 
small m whereas Fy(m) only increases for larger m. In 
all cases 70% of the dynamics is explained by the first ten 
modes, but only 35% of the dynamics taking place after 
the crack projects on the first three modes, whereas this 
fraction reaches 60% for the dynamics taking place be- 
fore or during the crack. A visual transcription of these 
numbers is provided by figure mH -fc). 
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FIG. 9: Projection of the dynamics restricted to the rigid 
subsets. Top: Fraction of the dynamics projected on the 
10 most significant mode: (F(10)) vs. r as defined in the 
text. The arrow on the right indicates the direction of in- 
creasing packing fraction (j} € [0.8417,0.8426,0.8440,0.8457]. 
Bottom: Average value of (_F(m)) vs. m for different 0. The 
error bars correspond to the standard deviations of F(m) and 
the dashed line is the prediction of Eq. ^ for a basis com- 
posed of random modes. 



These observations suggest that (i) the crack consid- 
ered here is really a micro-event in the sense that the 
dynamics after the crack still projects on a small amount 
of modes (here of the order of ten), (ii) at least for such 
a micro-crack, there is a selection of directions in phase 
space along which the cracks occurs. 



V. DISCUSSION AND CONCLUSION 

This paper is a first attempt to describe the dynamics 
of granular media in the jammed, glassy region in terms 
of "modes" , by applying a Principal Component Analysis 
(PCA) to the covariance matrix of the position of individ- 
ual grains. This is perfectly justified, and gives sensible 
results, in a regime of time/densities such that the aver- 
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FIG. 10: Analysis of a relaxation event within a sub-region, 
following a quiet period 4> = 0.8417. Above: (a) Cq{t) indi- 
cates a sudden collective motion of the particles, the so-called 
"crack". The four arrows enumerated from 1 to 4 indicate 
the bounds of the time intervals, where the real dynamics 
is considered. The quiet period during which the modes are 
computed is indicated by the straight horizontal (red) line 
just preceding the crack, (b) F(m) vs. m/2N for the three 
dynamics defined above; the vertical blue dotted line indi- 
cates m = 3. (c): On top of each other the real dynamics 
during the three periods of time and its recomposition with 
the 3 most significant modes 



age position of the particles is approximately constant, 
that is, varies less than the typical fluctuations them- 
selves, otherwise both the reference configuration and co- 
variance matrix itself evolve with time. The time scale 
over which the reference configuration can be considered 
as stable also depends on the system size, since in an in- 
finite system some rearrangement takes place somewhere 
in the system at each instant of time. 

For small enough times/system sizes, or at high enough 
packing fractions, this stability criterion is approximately 
fulfilled and the spectral properties of the covariancc 
matrix reveals large, collective fluctuation modes that 
cannot be explained by a Random Matrix benchmark 
where these correlations are discarded. The existence of 



these collective modes is expected from the results of [1^ 
that established the existence of dynamical correlations, 
which diverge as the system reaches its rigidity transition 
(f)j. The analysis in terms of eigenmodes provided here 
confirm that the slow, large scale dynamic structures ap- 
pear when — !■ , that explain a substantial fraction of 
the dynamics. 

We then attempted to find some link between the 
softest modes of the covariance matrix during a certain 
"quiet" time interval and the spatial structure of the re- 
arrangement event that ends this quiet period. In order 
to do so, we first tried to identify well-defined "cracks" 
that would lend themselves to such an analysis. This 
proves to be exquisitely difficult: the rearrangements are 
made of micro-cracks of all amplitudes, that span larger 
and larger regions of the system as cj) ^ and that 
are at the origin of the superdiffusive. Levy flight char- 
acter of the motion found in [l3|. In spite of this diffi- 
culty, we have succeeded in identifying some "rigid sub- 
sets" where well characterized cracks appear. The mo- 
tion during these cracks is indeed well explained by the 
soft modes of the dynamics before the crack. However, a 
more systematic analysis should be undertaken because 
we do not know at this stage whether the identification 
of these rigid subsets induces a strong selection bias on 
the nature of the cracks themselves. In the hypothesis 
where for a majority of cracks we could even not define 
the precursor modes then these soft modes would not be 
relevant to understand the dynamical evolution of the 
system. We believe that this is increasingly the case as 
one approaches the rigidity transition, where self-similar 
micro-cracks of all scales become overwhelming. In that 
eventuality, the analysis in terms of mode would only be 
useful to characterize the rigidity of amorphous granular 
systems, for dense enough packings where the rigid sub- 
sets remain dominant. In all cases, we believe that the 
methodology presented here will motivate and buttress 
further work in that direction. 
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